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REACTIVE  GAS  PHASE  COMPRESSION  DUE  TO  SHOCK- 
INDUCED  CAVITY  COLLAPSE  IN  ENERGETIC  MATERIALS 


Linhbao  Tran 

U.S.  Army  Research  Laboratory 
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Abstract.  A  mesoscale  simulation  is  carried  out  to  examine  shock-initiation  due  to  gas 
phase  reaction  at  site  of  cylindrical  pore  within  an  HMX  crystal.  The  focus  here  is  to 
investigate  viscoplastic  heating  with  gas  pore  compression  that  leads  to  chemical  reactions 
within  the  gas  phase.  Systems  of  conservation  laws  for  both  solid  and  gas  phases  are 
solved  along  with  species  conservation  from  a  reduced  set  of  chemical  kinetic  model. 
Mass,  momentum,  and  energy  transfer  between  phases  are  applied  explicitly  at  the  solid-gas 
interface  using  physical  boundary  conditions,  thus  avoiding  empiricism  of  mixture 
multiphase  formulation.  These  transfer  processes  are  critical  in  the  Mach  stem  formation 
around  the  collapsing  reacting  gas  pore. 


INTRODUCTION 

The  initiation  phenomenon  in  a  heterogeneous 
energetic  material,  such  as  a  plastic  bonded 
explosive  (PBX)  is  dominated  [1]  by  the  occurrence 
of  hot  spots  at  density  discontinuities  within  the 
material.  Initiation  of  an  energetic  material  can 
occur  when  an  impulse  delivered  to  a  high 
explosive  material  evolves  into  a  detonation  wave. 
Growth  to  detonation  depends  on  the  formation  of 
localized  regions  of  elevated  thermal  energy,  or  hot 
spots  that  have  temperatures  higher  than  the  bulk 
temperature  expected  from  shock  heating.  When 
sufficient  thermal  energy  is  generated  locally, 
ignition  occurs  and  subsequent  chemical  energy 
release,  having  overcome  energy  dissipation  in  the 
form  of  heat  conduction,  leads  to  progressive  shock 
strengthening  until  detonation  conditions  are 
achieved.  This  concept  was  first  proposed  by 
Bowden  and  Yoffe  [1], 

A  number  of  mechanisms  for  the  formation  of 
hot  spots  have  been  proposed  over  the  years,  such 
as  frictional  heating  between  adjacent  grains  [2-3], 


hydrodynamic  jetting  [4],  viscoplastic  void  collapse 
[5-7],  and  adiabatic  gas  pore  compression  [1,8-10]. 
The  applicability  of  these  mechanisms  depends  on 
the  strength  of  the  imposed  shock,  size  of 
microstructure  defects,  material  properties,  etc.  that 
initiates  the  explosion.  Regimes  of  applicability  can 
be  identified  roughly  by  the  criterion  proposed  by 
Khasainov  et  al.  [6]. 

In  the  current  study,  we  investigate  the  effect 
of  viscoplastic  heating  and  gas  pore  compression  of 
an  idealized  cylindrical  pore  due  to  the  passage  of  a 
shock.  The  formulation  involves  solving  both  solid 
and  gas  phases  separately  using  a  sharp  interface 
technique  and  shock  capturing  methodology.  The 
coupling  of  the  two  phases  is  carried  out  by 
applying  explicitly  the  physical  boundary  condition 
at  the  interface.  A  reduced  chemical  kinetic  model 
[11-12]  is  implemented  that  includes  three  reactions 
with  four  lumped  chemical  species. 

FORMULATION 

Governing  Equations 


The  two-dimensional  mass,  momentum,  and 
energy  conservation  equations  with  additional 
transport  equations  for  equivalent  plastic  strain  and 
deviatoric  stresses  in  the  solid  phase,  are  given 
generally  as: 

Qt+F(Q)x+G(Q)y=S{Q)  (1) 

where  vector  of  conserved  variables  is  Q , 
convective  flux  vectors  are  F(Q)  and  G(Q) ,  source 
term  is  S(Q) ,  and  subscripts  are  derivatives  with 


Where  N  is  the  number  of  species,  the  subscript 
i  is  for  i"'  species,  h,  is  the  species  enthalpy,  YI  is 
the  mass  fraction,  and  co ,  is  the  mass 
production/destruction  rate. 

In  the  above  source  terms  for  the  gas  phase  we 
assume  that  the  Lewis  number  is  identically  one, 

i.e.  yaD,  =  — — ,  where  D,  is  the  binary  diffusion 

C  p,i 

coefficient,  C  .  is  the  specific  heat,  and  /r  is  the 


respect  to  temporal  (t)  and  spatial  directions  (x  and 

y)- 

The  conservative  variables  and  convective 


thermal  conductivity.  Diffusion  velocities  are 
obtained  by  assuming  Fick’s  law  for  mass 
diffusion.  Here  we  neglect  the  Soret  and  pressure 


fluxes  for  gas  phase  are: 

Q  =  {p,  pu,  pv,  pE) 

(2) 

gradient  effects, 
U,=-yYx  and  Vj 

F(Q )  =  [pu,  pu 2  +  p,  puv,  u(pE  +  />)} 

(3) 

G(Q)  =  [pv,  puv,  pv 2  +  p,  v(pE  +  p)\ 

(4) 

In  the  solid  phase, : 
diffusion  velocities  are 

For  the  solid  phase,  the  deviatoric  stresses  are 
expressed  in  Jaumann  rate  [13]  form: 

Q  =  \p,  pu,  pv,  pE,  psp ,  psxx ,  ps)y ,  psxy } 

F(Q)  =  {pu,  pu 2  +  p,  piiv,  u(pE  +  p),  pus 


pus ^  pus pus ^ 


(5) 

(6) 


G(Q)  =  {pv,  puv,  pv~  +  p,  v(pE  +  p\pvs 


p’ 


(7) 


PVSxx’PVSyy’PVSxy, 

The  dilatational  component  of  the  stress  is 
obtained  from  the  equation  of  state  and  is  presented 
later.  Symbols  used  here  follow  standard 
convention. 

In  addition  to  the  usual  viscous  stresses  and 
thermal  conduction  in  the  source  terms,  there  are 
enthalpy  fluxes  due  to  species  diffusion  and  energy 
generation  due  to  chemical  reactions  in  the  energy 
equation.  For  the  gas  phase, 

0 

4  p 


s(e)= 


[us 


-PU,*x+—V,xy+PU,yy 

P  4 

PV,XX  +yM.V’  +^FV,yy 

+  vs  J ,x+(us*y  +vsyy\y  - 


N  \  N 

PlLKYyi 

i=i  7  v  /=i 


(8) 


—  Y 

Y  ’>y 


for  solid  phase  is  therefore  given  as: 

0 


S  +  V 

xx,  x  xy,y 


(9) 


s(e)= 


s  +s 

xy,x  _  yy,y 


\US  vv  +  VS 


),  +{“S*y  +VSyy),y  ~ 

+  T,yy  )-  X  uJh 
_  1=1 

jj yyi) 


3sxxI.  +  2nxysxy+2G\ 
3sxyl.  +  2nyxsxy  +  2G 


v  „  -S 


3s„S  +  ®  +  2G 


^iU,y+V,  J 


(10) 


Where <D=Q  s  +Q  s  —  Q  s  -Q  s  , 

xx  xy  xy  yy  xy  xx  yy  xy  ’ 

'Z  =  ^uii,D?  is  the  plastic  strain  rate  tensor,  and 


D.iy  =  -j  {uUj  +  u jj  )  is  the  spin  tensor. 


In  addition  to  the  above  system  of  equations, 
the  species  conservation  equations  are  solved  in 
each  phase,  given  as: 

{pYi),,+{pYiu),x+{pYiv),y  = 

{pDYlx)x+{pDYli)y+coi  <lh 

Mass  production  rates  are  determined  by 
phenomenological  chemical  kinetic  expressions. 


Assuming  that  the  reactions  are  thermodynamically 
simple: 


M 
/= 1 


KWJ  J 


(12) 


Where  M  is  the  number  of  reactions,  Wt  is  the 


molecular  weight  of  i"  species,  and  v'n  and  u",  are 

the  stoichiometric  coefficients  of  the  reactant  and 
product,  respectively. 

The  reaction  rate  is  based  on  an  Arrhenius 
form,  having  unit  of  s’1: 


kf  =  A  exp 


~Ea 

R..T 


(13) 


Where  A  is  the  pre-exponential  factor,  Ea  is  the 
activation  energy,  and  Ru  is  the  universal  gas 

constant.  The  expression  for  i,h  species  enthalpy 
is: 

hi  =  A/?;  ,-  +\lo  C p  idT  (14) 


Where  A h%  is  the  enthalpy  of  formation  for 

r  species  at  reference  temperature,  To ,  taken  as 
298  K. 


Chemical  Kinetic  Model 

A  reduced  kinetic  model  proposed  by  McGuire 
and  Tarver  [11]  is  implemented  to  model  the 
reactions  in  the  condensed  and  gas  phases.  We  are 
interested  only  in  how  qualitatively  distinct  species 
of  a  reduced  kinetics  model  manifest  themselves 
during  hot-spot  formation  and  the  subsequent 
ignition  of  the  reactive  material  in  a  viscoplastic 
void  collapse  process. 

The  thermal  decomposition  of  HMX  can  be 
modeled  as  a  three  stage  process  [11]: 

1  2  3 


HMX  — >  fragments  — >  intermediate  gases  — >  final  gases 


C4H8N808  CH2NN02 

Or  equivalently, 
A  ->  B  -> 


CH20,  n2o, 
HCN,  HN02 

2C  ->  D 


N,,  H20, 
C02,  CO 


The  three  major  groups  of  irreversible  rate¬ 
controlling  steps  in  the  HMX  decomposition 
process  are  (1)  unimolecular  endothermic  breaking 
of  C-N  and  /  or  N-N  bonds  to  produce  methylene 
nitramine  and  other  higher  molecular  weight 
fragments;  (2)  weakly  exothermic  unimolecular 
decomposition  of  these  condensed  phase  fragments 
to  produce  intermediate  gases:  CH20  +  N20  or 


HCN  +  HN02;  and  (3)  highly  exothermic 
bimolecular  gas  phase  reactions  between  these 
intermediate  gases  to  produce  stable  reaction 
products  N2,  H20,  C02,  CO,  etc.  Note  that  of  these 
three  mechanisms  the  first  occurs  in  the  condensed 
phase,  the  last  two  reactions  occur  in  gas  phase. 
For  vapor  reactant  (species  B)  to  exist  in  the  gas 
pore,  formation  of  species  B  in  the  solid  phase  is 
allowed  to  advect  into  the  gas  phase  through  the 
boundary  condition  apply  at  the  gas  phase’s  ghost 
cells. 


Constitutive  Relations 

The  equations  governing  the  material 
deformation  appropriate  for  high  strain  rate 
applications  can  be  formulated  by  assuming  that  the 
volumetric  or  dilatational  response  is  governed  by 
an  equation  of  state  while  the  shear  or  deviatoric 
response  obeys  a  conventional  flow  theory  of 
plasticity.  The  equation  of  state  for  solid  phase  is 
of  Mie-Griineisen  form: 


p(e,v)=r^+f(v) 

The  Griineisen  parameter  is  defined  as: 

■TqPo 

V  p 


(15) 


r.nfX 

de 


Poco  V 


1  1 


V  V, 


L ,  V  =  Yp ,  and 

if 

V<V0 

2V  J 

] 

[f 

v  >V0 

0  J 


(16) 

The  yield  strength  is  assumed  to  be  temperature 
dependent. 


,(r)=. 


T,„-T 
~T0  J 


Where  the  melting  temperature  [13]  is: 

'  2(r.  -  xl 


T,Xp)=Tm, 


1+- 


p 


\p~Po ) 


V  J 

and  solid  phase  viscosity  is: 

[31.0  Pa-s  if  T  <Tm 

|3.1xl0-7exp(7800/)pa.s  ifT>Tm 


P  =  ‘ 


(17) 


(18) 


(19) 


Parameters  for  HMX  material  properties  taken 
from  [14]  are  given  to  be:  p0  =  1900  kg/m3,  C  = 
1031  J/kg-K,  Tmo  =  520  K,  s  =  2.38,  oyo  =  0.37 


GPa,  G  =  10.0  GPa,  T0  =  1.1,  and  c0  =  2650  m/s. 
For  chemical  and  thermal  properties,  refer  to 
Tarver  et  al.  [12]. 

For  gas  phase,  although  the  JWL  equation  of 
state  is  not  accurate  at  low  pressure,  we  opt  to  use 
this  form  as  we  expect  high  gas  pressure  at  the  later 
stage  of  the  pore  collapse  process,  more  critical  for 
the  problem  of  interest  here.  The  JWL  equation  is 
given  as: 


p(e,V)=A 


1- 


CO 

Ry 


exi 


(20) 


B\ 


1 - — 

v 


exp  (~R2V)+^r 


V 


Here  V  =  p0j  p  and  A  =  38.0651  Mbar,  B 


1.2948  Mbar,  R,  =  7.70,  R2  =  2.40,  co  =  0.33.  The 
specific  heat  for  gas  phase  is  based  on  mass 

N 

fraction  averaged:  Cp  =  X  Yt  C pi  ,  and  similarly  for 

1=1 

N 

thermal  conductivity:  k  =  £  YjKi 

i=l 


where  the  interface  normal  velocity  is  vI  n  and  m"  is 

the  interfacial  mass  flux. 

The  balance  of  momentum  must  be  ensured 
across  a  solid-gas  interface  which  has  zero 
thickness.  This  requires  that  the  sum  of  traction  for 
the  two  materials  vanish: 

Cgas  +  "t solid  =  0  (22) 

For  the  gas  phase,  this  traction  tensor  reduces 
to  just  the  hydrostatic  pressure  and  viscous  stresses. 
Kinematics  boundary  conditions  are  also  applied 
for  velocities  across  the  interface. 

Temperature  continuity  at  the  interface 
requires: 

T,  =  T.,  =  TgJ  (23) 

with  interface  temperature,  7j  obtained  from  the 
balance  of  heat  fluxes  due  to  conduction  of  the  two 
phases.  Only  reactant  species  B  is  allowed  to  be 
advected  into  the  pore  in  form  of  reactant  vapor. 
The  change  in  internal  energy  for  gas  phase  due  to 
incoming  species  B  is:  YBh°f  B 


Initial  Conditions 

The  initial  conditions  describe  an  equilibrium 
state  between  the  gas  phase  and  condensed  phase. 
The  pore  is  initially  considered  to  be  completely 
tilled  with  final  product  gas,  species  D. 
Temperatures  are  given  as: 

T  =  T  =T 

e  g  o 

in  the  condensed  phase,  species  mass  fraction  are 
given  to  be: 

ya  =1.0;  yb=yc  =  yd=  o.o 

And  in  the  gas  phase: 

Yd  =1.0;  Ya=Yb  =  Yc=  0.0 
Pressures  and  densities  are  taken  to  be: 

Pg  ~  P  o  ’  P  g  ~  Pg,o 

where  the  subscript  “o”  denotes  the  reference  state 
with  values: 

Ta  =  298 K,pa  =1.0 Bar,  and pa g  =  0.897 kg/m3 

Interfacial  Boundary  Conditions 

The  gas-solid  interface  boundary  conditions 
can  be  derived  from  conservation  laws  by 
considering  a  control  volume  surrounding  a 
segment  of  the  interface.  The  jump  condition  in 
mass  conservation  is  given  as: 

™  ”  =  Ps  (V  V,  "  VU,  )  =  Pg  (Vg,„  -  Vi,n  )  (2 1 ) 


Numerical  Method 

Spatial  discretization  with  a  second  order 
convex  ENO  scheme  [15]  and  a  third  order  Runge- 
Kutta  for  time  stepping  were  implemented.  Near 
the  material  interface,  fluxes  are  modified  to 
account  for  the  interface  boundary  conditions  [16]. 
The  presence  of  stiff  chemical  source  terms 
requires  the  usage  of  an  operator  splitting  technique 
[17].  The  resulting  ODEs  are  nonlinear  and 
discretized  with  a  fully  implicit  first  order  Euler 
scheme.  For  accuracy  and  stability  of  the  method, 
the  time  step  for  the  governing  equations  is 
adaptive,  with  a  lower  limit  for  CFL  value 
according  to  the  degree  of  stiffness  in  the  source 
term.  Subcycling  for  the  ODEs  is  also  employed 
for  enhanced  computational  speed. 

In  the  solid  phase,  viscoplastic  regularization 
algorithm  [18]  is  applied  to  ensure  the  correct 
admissible  stress  states. 

The  interface  is  captured  using  level-set 
methodology.  The  transport  for  the  level-set 
equation  is  solved  using  3rd-order  ENO  scheme 
along  with  3rd-order  Runge-Kutta  time  stepping. 

RESULTS 

The  setup  for  the  problem  of  interest  here  is  a 
cylindrical  void  5  pm  in  diameter,  centered  at 


(0,12.5)  jum  in  a  (12.5  jumx 25 /m)  domain. 
Symmetry  boundary  conditions  are  applied  on  left 
and  right,  with  inflow  and  outflow  boundary 
conditions  on  the  bottom  and  top  boundaries. 
Initially  a  piston  is  driven  into  the  HMX  crystal 
from  the  bottom  by  applying  a  material  velocity  of 
1  km  Is,  for  a  pulse  of  1  ns ,  see  insert  sketch  of 
Fig.  1. 

Pressure  profiles  of  three  probes  positioned 
near  the  right  boundary  are  shown  in  Fig.  1.  The 
pore  would  experience  an  average  shock  pressure 
of  about  6  GPa  with  a  rise  in  temperature  to  around 
700  K  due  to  shock  heating.  The  lower  surface 
material  velocity  is  calculated  to  be  around 
1.7  km  I  s  .  For  a  hydrodynamic  pore  collapse 
process,  the  characteristic  time,  Th  =  d  /  2 U  m  is 
about  1.5  ns .  In  our  calculation,  the  collapse  time 
is  around  2.2  ns  . 


FIGURE  1.  Solid  phase  pressure  profiles. 

As  the  shock  front  propagates  through  the  pore, 
there  is  a  rapid  decrease  in  peak  shock  pressure 
from  about  6.3  GPa  to  5.4  GPa,  within  a  travel 
distance  of  only  5  jam  .  There  is  also  a  clear 
change  in  slope  toward  the  top  part  of  pressure 
profiles.  The  first  effect  is  due  to  rarefaction  waves 
from  behind  attenuating  the  shock  pressure.  The 
second  effect,  the  change  in  pressure  slope,  is  due 
to  viscoplastic  regularization  formulation,  where 
stresses  are  allowed  to  relax  to  the  yield  surface  and 
the  rate  depends  on  material  viscosity.  Note  that 
this  is  not  a  two  wave  structure  (an  elastic  precursor 
followed  by  the  plastic  wave,)  as  the  shock  wave  is 
strong  enough  to  overrun  the  elastic  wave. 

Grid  convergence  in  the  solid  phase  is 
demonstrated  (Fig.  2)  utilizing  four  different  mesh 
densities  (50x100,  100x200,  150x300,  and 


200x400).  Pressure  profiles  at  three  locations  (Fig. 
1)  in  the  domain  show  only  a  slight  difference  for 
the  two  fine  meshes.  We,  therefore,  assume  the 
finest  resolution  at  200x400  as  our  grid  independent 
solution. 


FIGURE  2.  Grid  independent  study  -  Solid 
phase  pressure  profiles. 

When  we  plot  the  profiles  for  maximum  gas 
phase  temperatures  (Fig.  3),  the  solution  does  not 
yet  seem  to  converged.  For  early  part  of  the 
profiles,  up  to  about  2.5  ns ,  the  solutions 
demonstrate  convergence,  however  at  later  time, 
temperature  of  the  finer  mesh  begins  to  rise  earlier 
than  that  of  the  coarse  meshes.  From  2.5  ns  to 
about  3.5  ns,  although  the  differences  are 
noticeable,  they  are  not  too  drastic.  To  this  point, 
much  of  the  temperature  rise  is  due  to  gas 
compression  within  the  pore  as  well  as  mass 
transfer  into  the  pore  in  form  of  vapor  species  B. 
Note  that  there  exists  gas  shocks  present  within  the 
pore. 


FIGURE  3.  Grid  independent  study  -  Maximum 
gas  phase  temperature. 


Above  3.5  ns,  the  solutions  from  the  four 
meshes  show  drastic  differences  in  term  of  the 
steepness  of  temperature  rise  and  peak  magnitude. 
This  is  due  to  the  fact  that  we  have  a  chemically 
reactive  system  that  begins  to  react  given  sufficient 
rise  in  temperature. 

For  a  truly  converged  solution,  one  would  need 
to  have  a  time  step  that  at  least  has  the  same  order 
of  magnitude  of  the  time  step  of  the  fastest  reaction 
within  the  system.  This  would  require  orders  of 
magnitude  higher  resolution  and  is  clearly  too 
computationally  expensive  at  present. 

A  similar  conclusion  is  reached  when  we  plot 
average  gas  pressure  versus  time  (Fig.  4).  For  the 
coarsest  resolution  (50x100),  the  rise  in  pressure  is 
only  an  order  of  magnitude  above  the  initial 
condition,  whereas  for  the  finest  mesh,  pressure  rise 
can  reach  up  to  1000  bars.  The  drop-off  toward  the 
end  is  due  to  the  fact  that  there  is  not  enough 
resolution  in  the  collapsing  pore. 

Even  with  this  large  rise  in  gas  pressure,  it  is 
still  lower  than  the  yield  strength  of  the  solid  phase 
to  cause  a  reversal  of  plastic  flow.  Flowever,  as 
will  be  shown  later,  this  provides  enough  resistance 
to  slow  down  the  pore  collapse  process  for  the 
critical  gas  phase  energy  release  to  transfer  to  the 
solid  phase.  For  subsequent  discussions,  we  will 
use  the  results  from  the  200x400  mesh  size 
calculations. 


FIGURE  4.  Grid  independent  study  -  Average 
gas  phase  pressure. 

Figure  5  shows  profiles  of  maximum  gas  and 
solid  phase  temperature,  along  with  average  mass 
fraction  of  gaseous  species  within  the  pore.  The 
initial  rise  in  solid  phase  temperature  up  to  about 
1.2  ns  is  due  to  shock  heating.  Thereafter,  until 
about  2.5  ns  ,  the  drop  off  is  due  to  the  finite  shock 


pulse  width  giving  rise  to  shock  attenuation.  Gas 
pore  temperature  starts  to  increase  when  the  shock 
hits  the  lower  surface  of  the  pore,  around  2.1  ns  ,  up 
to  about  3.6  ns.  The  increase  is  also  due  to 
interfacial  mass  flux  coming  from  the  solid  phase. 

Note  that  from  about  2.8  ns  to  around  4  ns , 
the  maximum  gas  phase  temperature  is  higher  than 
that  of  solid  phase.  The  higher  rise  in  temperature 
is  due  mostly  to  gas  phase  compression  along  with 
interface  mass  flux  coming  from  the  solid  phase. 
There  is  negligible  exothermic  reaction  at  this 
temperature. 


FIGURE  5.  Maximum  temperature  in  solid  and 
gas  phases  in  lower  profiles.  Average  gaseous 
species  mass  fraction  in  the  upper  profiles. 

Just  beyond  4  ns ,  there  is  a  rapid  rise  in  both 
gas  and  solid  phase  temperature  with  gas  phase 
temperature  reaching  4000  K  around  4.2  ns  .  The 
solid  phase  temperature  rises  above  that  of  the  gas 
phase,  due  to  the  formation  of  a  solid  phase  shock 
stem  near  the  surface,  a  region  where  solid  phase 
has  fully  undergone  conversion  from  species  A  to 
B,  i.e.,  does  not  require  any  more  energy  for 
endothermic  reaction.  The  third  reaction  starts  to 
kick  in  around  4.2  ns  as  seen  in  the  abrupt  rise  of 
species  D  mass  fraction,  along  with  the  abrupt  rise 
of  gas  temperature  (maximum  around  5500  K). 
One  can  also  see  the  rapid  consumption  of  species 
B  and  C  as  well.  From  this  point  on,  the  rate  of 
vapor  reactant  supplied  to  the  gas  pore  cannot  keep 
up  with  the  consumption  rate.  Therefore,  for  gas 
phase  chemical  energy  release,  the  rate  limited 
process  is  the  vapor  reactant  mass  flux  at  the  gas- 
solid  interface,  as  expected. 

Pore  collapse  under  the  current  dynamic  shock 
loading  condition  is  not  a  symmetric  process  about 
the  center  (Fig.  6).  This  holds  especially  true  if  one 


considers  the  fact  that  gas  present  within  the  pore 
does  not  immediately  reach  an  equilibrium  state  for 
the  duration  of  pore  collapse  process.  As  the 
planar  shock  impacts  the  lower  pore  surface,  the 
interface  begins  to  accelerate,  but  at  the  same  time, 
the  gas  is  being  compressed  around  this  region  as 
well.  Although  gas  pressure  is  still  a  few  orders  of 
magnitude  smaller  than  that  of  solid  pressure,  it 
provides  a  resistance  to  somewhat  slow  down  the 


FIGURE  6.  Close-up  view  of  gas  phase  pressure 
and  velocity  field  show  shock  formation  within 
the  pore. 

Once  the  shock  in  the  solid  phase  moves  up  to 
the  side  of  the  pore,  the  transverse  impingement  on 
the  side  surface  causes  a  rapid  collapse  process, 
giving  substantial  rise  in  temperature  due  to 
viscoplastic  heating  compared  to  bulk  shock 
heating.  This  is  critical  as  the  temperature  in  this 
region  can  be  high  enough  for  solid  phase 
transformation,  the  first  endothermic  reaction  to 
take  place  which  will  supply  the  pore  with  reactant 
vapor.  Another  interesting  feature  resulting  from 
irregular  pore  collapse  is  shear  banding  mechanism 
as  also  shown  in  the  figure.  This  does  give  rise  to 
localized  thermal  field  slightly  away  from  the  pore 
surface. 

Contour  plots  of  species  B  (left)  and  C  (right) 
mass  fractions  (Fig.  7)  show  the  depletion  of 
species  B  and  the  production  of  species  C.  One 
can  see  clearly  the  complete  destruction  of  species 
B  near  the  solid-gas  interface  and  the  maximum 
production  of  species  C  in  the  same  region.  This  is 
in  essence  a  flame  front  due  to  the  exothermic 
reaction.  Elsewhere  in  the  pore,  the  temperature  is 
not  high  enough  for  spontaneous  2nd  reaction  to 
take  place,  although  there  is  a  continuous 


production  of  species  C.  Production  of  species  D, 
the  last  exothermic  reaction,  with  locally  much 
smaller  reaction  rate,  is  minimal  at  this  point  in 
time.  In  the  region  around  the  interface  where  gas 
phase  reactions  occur,  the  energy  is  being 
continuously  transferred  to  the  solid  phase.  This 
results  in  a  complete  conversion  of  species  A  to 
species  B  as  shown  in  the  left  of  figure  7. 


FIGURE  7.  Close-up  view  shows  destruction  of 
species  B  (left)  inside  the  pore  as  well  as 
production  of  species  C  (right)  -  Notice  thin 
layer  adjacent  to  pore  surface. 

In  figure  8,  we  show  a  time  sequence  of 
pressure  contours.  From  the  shape  of  the  gas  pore, 
the  collapse  process  is  initially  not  symmetric  about 
the  pore  center,  and  in  later  time,  when  the  gas  pore 
pressure  equilibrated,  it  becomes  symmetrical.  As 
the  shock  plane  impacts  the  lower  pore  surface,  it 
produces  an  essentially  free  surface  since  the  initial 
gas  pressure  is  negligible  compare  to  shock 
pressure.  The  rarefaction  wave  created  by  this  free 
surface  travels  back  into  the  compressed  material, 
attenuating  the  shock  profile  close  to  the  pore 
surface.  In  the  middle  plot  on  second  row,  a  Mach 
stem  (next  to  pore  surface)  starts  to  form  as  the  gas 
pore  begins  to  react.  This  is  close  to  the  location 
where  shear  banding  occurs.  However,  for 
calculation  without  chemistry,  shear  banding  alone 
is  not  a  significant  effect.  So  there  is  a  coupling  of 
the  release  of  chemical  energy  inside  the  pore  to  the 
solid  phase  energy  at  the  interface.  In  subsequent 
plots,  the  Mach  stem  grows  stronger,  becomes 
circular,  and  begins  to  detach  from  the  pore  surface. 
From  this  point  on,  the  Mach  stem  grows  stronger 
than  the  leading  shock  and  it  will  eventually  catch 
up  and  over-take  the  leading  shock.  Again,  the 
speed  of  the  shock  is  enhanced  due  to  the  medium 
already  being  compressed. 

Figure  9  shows  a  close-up  view  of  the  shock 


-  flame  front 


v  gas-solid  interface 


FIGURE  8.  Sequence  of  pressure  contours  show  Mach  stem  formation  due  to  gas  phase  chemical 
energy  release  around  the  collapsing  pore. 

detaching  from  the  pore  surface,  as  well  as  the  Mach  stem  is  of  order  10  GPa  ,  causing  the  pore  to 

velocity  vector  field.  The  pressure  field  behind  the  continue  to  collapse  since  the  gas  pressure, 


(although  3  orders  of  magnitude  higher  than  initial 
value)  is  not  enough  to  cause  a  reversal  in  plastic 
flow.  The  pore  shape  at  this  stage  is  rather 
symmetric.  Note  that  within  this  circular  shock, 
full  conversion  from  solid  species  A  to  B  is 
completed  (Fig.  7).  This  in  turn  continuously 
supplies  the  pore  with  fresh  vapor  reactant  in  form 
of  species  B. 


FIGURE  9.  Close-up  view  of  shock  formation  as 
it  detaches  from  the  pore  surface. 

CONCLUSION  AND  REMARK 

Conclusion 

A  mesoscale  simulation  is  carried  out  to 
examine  shock-initiation  due  to  gas  phase  reaction 
at  site  of  cylindrical  pore  within  an  HMX  crystal. 
Systems  of  conservation  laws  for  both  solid  and 
gas  phases  are  solved  along  with  species 
conservation  from  a  reduced  set  of  chemical  kinetic 
model.  Mass,  momentum,  and  energy  transfer 
between  phases  are  applied  explicitly  at  the  solid- 
gas  interface  using  physical  boundary  conditions. 
It  is  found  that  these  processes  provide  the 
necessary  and  critical  bridge  for  chemical  energy 
generated  within  the  gas  pore  to  transfer  to  the 
solid  phase. 

Although  we  observe  a  flame  front  near  a 
region  of  the  solid-gas  interface,  the  interface 
regression  rate  is  limited  to  a  mass  diffusion 
velocity.  Therefore,  the  outward  propagation  of 
the  interface  due  to  surface  burning  is  small 
compare  to  the  pore  collapse  rate.  The  transfer  of 
thermal  energy  to  the  solid  layer  adjacent  to  the 
pore  gives  rise  to  temperature  which  is  sufficient 
for  the  continuing  first  endothermic  reaction  in  the 
solid  phase  to  take  place.  This  in  turn,  provides 
fresh  vapor  reactant  to  the  gas  pore,  necessary  for 
continuing  combustion.  The  feed  back  cycle  is 
limited  to  the  mass  flow  rate  of  vapor  reactant 


coming  into  the  pore.  The  gas  phase  chemical 
energy  release  transferred  to  the  solid  phase,  causes 
a  rise  in  the  solid  phase  pressure.  Mach  stem 
formation  is  observed  around  the  shear  banding 
region.  Although  shear  banding  is  not  the  cause  of 
shock  formation,  it  contributes  to  the  localized 
elevated  thermal  energy.  When  the  shock  gathers 
enough  strength,  it  detaches  from  the  pore  surface 
and  propagates  outward. 

With  respect  to  the  reversal  of  plastic  flow 
around  the  pore  observed  in  other  idealized 
spherical  viscoplastic  models,  this  study  shows  that 
the  gas  pore  pressure  never  reaches  values  above 
the  pressure  field  in  the  solid  phase  around  the 
pore.  Although  the  gas  phase  pressure  could  reach 
up  to  3  orders  of  magnitude  higher  than  initial 
value,  it  never  exceeds  the  surrounding  solid 
pressure.  In  theory,  the  gas  pressure  could  be  very 
high  as  the  pore  collapses,  however,  the  resolution 
of  the  pore  is  poor  towards  the  end  of  its  life  and 
therefore  numerical  dissipation  could  be 
significant.  On  the  other  hand,  the  physical 
reasoning  is  that  with  continuous  transfer  of  gas 
phase  chemical  energy  to  solid  phase,  the  solid 
phase  pressure  around  the  interface  also  rises, 
which,  it  appears  the  gas  pressure  could  not 
overcome. 

It  is  shown  in  this  study  that  viscoplastic  pore 
collapse  is  a  non-equilibrium  process.  Although 
the  characteristic  collapse  time  is  close  to  the 
hydrodynamic  mode  as  categorized  by  other  works, 
the  process  itself  is  not  hydrodynamic.  Previous 
works  rely  on  an  idealized  spherical  pore  to 
categorize  pore  collapse  regimes  with  gas  phase 
pressure  and  temperature  taken  from  bulk  values, 
similarly  for  the  solid  phase.  The  collapsing 
process  is  quite  complex  however  and  non¬ 
equilibrium,  competing  processes  must  be 
accounted  for. 

Remark  and  Future  Work 

Interface  energy  transfer  is  formulated  based 
only  on  the  balance  of  conduction  heat  fluxes 
between  the  gas  and  solid  phases.  It  is  shown  that 
the  gas  phase  temperature  is  well  above  1000  K  in 
the  later  stage  of  pore  collapse  process  when  the 
gas  shock  propagates  away  from  the  pore  surface. 
Radiative  heat  transfer  could  play  a  role  in 
increasing  the  solid  phase  temperature. 

Interfacial  mass  flux  which  accounts  for  the 
vapor  reactant  coming  into  the  pore  is  based  on 


convective  flux  applied  at  the  interface.  The 
available  vapor  reactant  is  taken  from  solid  phase 
and  applied  to  the  ghost  cells  of  gas  phase.  More 
rigorous  interfacial  mass  flux  treatment  is  needed. 

The  validity  of  the  JWL  equation  of  state  for 
the  initial  gas  phase  is  questionable.  The  initial  gas 
phase  pressure  buildup  although  not  critical,  is 
important  in  determining  how  the  pore  shape 
evolves.  This  plays  a  key  role  in  where 
viscoplastic  heating  is  dominant.  For  the  solid 
phase,  the  validity  of  the  Mie-Griineisen  equation 
of  state  for  a  re-shocked  material  is  unknown.  This 
issue  is  relevant  to  the  formation  of  the  Mach  stem 
in  the  shocked  solid  phase. 

In  term  of  numerics,  subcycling  of  the  ODEs 
along  with  adaptive  time  stepping  for  governing 
equations,  many  time  iterations,  0(10000)  are 
required.  This  could  give  rise  to  significant  round¬ 
off  error. 
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